%% Convolution;  The Ising Model
%

%% Convolution
% The convolution of two functions is defined as
%
%%
% $$f(x)\ast g(x)=\int_{-\infty}^\infty f(x)g(x-x')dx' $$
%
% The convolution theorem states that we can perform the convolution using
% Fourier transforms:
%
% $$ f \ast g = F^{-1}\{ F\{f\}*F\{g\}\}$$
%
% Where $F$ is the Fourier transform.
%%
% <https://en.wikipedia.org/wiki/Convolution_theorem>

% Lets try something basic

N=100;
X=linspace(0,10,N);
F = zeros(size(X));
F(1)=1;
F(round(N/2))=1;
G=exppdf(X);
figure
plot(X,F)
hold on
plot(X,G)
xlabel('X')
ylabel('Signal')
legend('F','G')

FoG=ifft(fft(F).*fft(G));
plot(X,FoG)









%% 1D Ising Model
% The Ising model considers a lattice of paricles that have up or down
% spin. The total energy is given by the interaction of nearest neighbors
% and interaction of the spin with an external field.
%
% $$H(\sigma)=-\sum_{<i,j>}J_{ij}\sigma_i\sigma_j-\mu\sum_jh_j\sigma_j $$
%
% Where the $<i,j>$ means adjacent particles, $h$ is
% the field, $\sigma$ is the spin, which can be either -1 or 1, and $J$ is
% the interaction energy.

%%
% <https://en.wikipedia.org/wiki/Ising_model>
clc
close all
clear all

Beta=1;         %1/KbT
J=0;            %Interaction Strength
B=2;            %Magnetic Field
N=40;           %Number of spins

%Start with random orientation
Spins=double(randn(N,1)>0);
Spins(Spins==0)=-1;

%Plot Spins
X=(1:N)';
Y=Spins==-1;
FigH=figure;
FigH.Position=[200,300,900,400];
quiver(X,Y,X*0,Spins,0,'linewidth',2)
axis([0 N+1 -1 2])
s=sprintf('Beta=%g, J=%g, B=%g',Beta,J,B);
title(s)

%Setup so we stop when figure is closed
FigH.CloseRequestFcn='Go=0;';
Go=1;

%Run Markov Chain Monte Carlo
NChain=1000000;
nn=0;
while Go
    nn=nn+1;    %counter for display
    
    %Make test jump
    
    %Pick site and flip spin
    ID=ceil(N*rand);
    SpinsTest=Spins;
    SpinsTest(ID)=-Spins(ID);
    
    %Energy from nearest neighbor interactions
    E_Spin=J*(Spins(1:end-1).*Spins(2:end)+Spins(1)*Spins(end));
    E_SpinTest=J*(SpinsTest(1:end-1).*SpinsTest(2:end)+SpinsTest(1)*SpinsTest(end));
    
    %Energy from Interactions with feild
    E_Field=-B*sum(Spins);
    E_FieldTest=-B*sum(SpinsTest);
    
    %Total Energy (Hamiltonian)
    H=E_Spin+E_Field;
    HTest=E_SpinTest+E_FieldTest;
    
    %Test acceptance probability for Metropolis-Hastings
    Alpha=exp(-Beta*(HTest-H));
    
    %Metropolis-Hastings rule
    if rand()< (Alpha)
        Spins(ID)=-Spins(ID);
    end
    
    ShowIter=1;
    %Show every ShowIter interations
    if nn/ShowIter==round(nn/ShowIter)
        figure(FigH);
        hold off
        Y=Spins==-1;
        quiver(X,Y,X*0,Spins,0,'linewidth',2)
        axis([-1 N+1 -3 4])
        title(s)
        hold on
        quiver(0,0,0,B,0,'linewidth',2)
        legend('Spins','Field')
        pause(.1)
    end
end
delete(FigH)
figure;
hold off
Y=Spins==-1;
quiver(X,Y,X*0,Spins,0,'linewidth',2)
axis([-1 N+1 -3 4])
title(s)
hold on
quiver(0,0,0,B,0,'linewidth',2)
legend('Spins','Field')
